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Abstract 

We report the first application of complex symmetric 
wavelets to the numerical simulation of a nonlinear signal 
propagation model. This model is the so-called nonlinear 
Schrodinger equation that describes, for instance, the evo- 
lution of the electric field amplitude in nonlinear optical 
fibers. We propose and study a new way to implement 
a global space-time adaptive grid, based on interpolation 
properties of higher-order scaling functions. 

1. Introduction 

The idea of using wavelets to perform numerical simu- 
lations of partial differential equations is not new (see, for 
instance, [1, 2] and references therein). The motivations 
come from the fact that wavelets provide a mathemati- 
cal representation that can resolve numerical difficulties 
due to singular phenomena. More exactly, wavelet prop- 
erties such as orthogonality and compact support, as well 
as exact representation of polynomials of a fixed degree 
in wavelet bases, allow efficient and stable calculation of 
regions with transient phenomena or strong oscillations. 
Up to now, wavelets that have been used for such stud- 
ies are the "classical" ones (real Daubechies wavelets [3], 
splines, Shannon and Meyer wavelets, etc.). 

In the present work, we are interested in using com- 
plex symmetric wavelets [4] as an interpolating tool in the 
numerical sampling process. To this aim, we will consider 
the nonlinear signal propagation model 

iu t + ^u xx + \u\u\ 2 = 0; AGlR*, (1.1) 

that is known as the nonlinear Schrodinger (NLS) equa- 
tion. Equation (1.1) describes, for instance, the evolution 
of the electric field amplitude u(t,x) in nonlinear optical 
fibers (in this particular application, t and x are space 
and time variables respectively). In this paper, we use 
the notation NLS"+" for A > and NLS"-" for A < 0. 
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The choice of the model (1.1) is interesting because it 
can lead to strong gradient phenomena. For the NLS"—" 
case, for instance, smooth localized initial conditions can 
evolve toward a typical breaking wave phenomena that ex- 
hibits very high frequency oscillations [5, 6]. In addition, 
breather solutions (also known as the bound N-soliton so- 
lutions in optics) of the NLS"+" equation, shows periodic 
peaking of the signal amplitude (see, for instance, [7] and 
references therein) and can be used to test the ability of 
an adaptive algorithm to increase and decrease the grid. 

The article is organized as follows. In Section 2, we 
recall the discrete wavelet theory through the multires- 
olution analysis and give an example of complex sym- 
metric wavelets. In Section 3, we study a simple way to 
implement a global space-time adaptive discretisation in 
the split-step Fourier method based on the recomposition 
properties of higher-order scaling functions. This globally 
adaptive algorithm will be tested on typical simulations 
(optical breaking waves and N-soliton solutions). 

2. Complex symmetric wavelets 

Wavelets can easily be introduced using the multires- 
olution analysis [8]. 

The basic equation of the multiresolution theory is 
the scaling equation that establishes a link between the 
underlying symmetries of the wavelet theory: dilations 
and translations. Given a set of complex coefficients 
ak,k £K, the scaling equation 

ip(x) = 2^2a k ip(2x - k), x G IR (2.1) 

k 

with the normalization J tp[x) dx = au = 1, define a 
scaling function <p(x). 

The multiresolution analysis is a particular decompo- 
sition of the space of square integrable functions i 2 (IR) 
into a chain of closed subspaces . . . Vj—i C V 3 C Vj+i C 

. . ., j G 7L, where the set {<pj,k(x) = 2? ip(2 3 x — k), k eTL) 
is a basis for Vj . The index j represents the scale to 
which the basic scaling function <p(x) is dilated. For a 
given scale, the space Vj is generated by the translations 
of ip h o(x) by the amount k2~ 3 . 
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Multiresolution also allows to decompose i 2 (IR) as 



i 2 (lR) 



v No e w j 



(2.2) 



3>N 



where Wj is the orthogonal complement of Vj in Vj+i 
(Vj + i = Vj ® Wj) and No is the coarsest resolution level. 

For a given scale j, Wj is generated by the set 
{ipj,k(x) = 2iip(2 J x -k),k £ 7L}, where the analysing 
wavelet ip( x ) satisfies 



■ip(x) 
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k), x £ 1R, 



(2.3) 



with bk = ( — l) k ai-k. 



Therefore, any function f(x) of 
i 2 (]R) can be expanded as a linear combination of trans- 
lates of the scaling equation at the scale No and the 
translates of the associated wavelet at the smaller scales 
No + 1,N + 2, as 
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w J,k1pj,k(x), (2.4) 



where w h k are the wavelet coefficients (that are non- 
negligible only in high gradient regions). In practice, ex- 
pansion (2.6) is computed by fixing a maximal resolution 
N for which the wavelet coefficients WN,k are negligible, 
that is, 

f(x) = ^ VN,k<PN,k(x). (2.5) 
k 

In this work, we will consider Daubechies' wavelets 
and, in particular, the recently investigated complex sym- 
metric cases [4]. Important numerical features of wavelets 
(compact support, orthogonality and regularity) can be 
implemented in the multiresolution analysis as the fol- 
lowing. Compact support requires a finite number of 
scaling coefficients a* in (2.1). Orthonormality condi- 
tion imposes this number to be even. In this work we 
choose k = — J, — J + I, J + I , which makes <p(x) and 
ip(x) non-zero inside the interval [—J, J + I]. Symme- 
try of the scaling function <p(x) with respect to x = 1/2 
is possible only for complex scaling coefficients. Due to 
the localisation properties of ipN,k(x) around the sam- 
pling points XN,k = 2^+1 ' Taylor expansion of fix) in 
VN,k = f <pN,k(x)f(x)dx then leads to 



VN,k 
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-N/2 
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(2.6) 



Coefficients v and w at lower resolutions are obtained from 
equations (2.1) and (2.3). Finally, regularity amounts to 
require that the first polynomial terms of the Taylor ex- 
pansion of a function / are mapped exactly in Vj . For a 
given order J, we can require such condition up to poly- 
nomials of degree J. An example of complex symmetric 
scaling functions and wavelets that result from the above 
requirements is given in Figure I for J = 4. 

In the following section, we will use this wavelet to 
detect strong gradients in the signal u(t, x) and to dy- 
namically change the sampling rate. 
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Figure 1: Real and imaginary parts of <p{x) and tp(x) 
for the complex symmetric case J = 4 



3. A globally adaptive algorithm 

One of the most popular numerical scheme to solve 
NLS-type equations is the so-called "symmetrized split- 
step Fourier method" (see, for instance, [7] and references 
therein). This is essentially a pseudo-spectral algorithm 
that splits the linear and nonlinear effects over one time 
step. Because nonlinear terms are more easily calculated 
in the physical space, only the linear part is decomposed 
in its Fourier components using a FFT algorithm. 

We now want to demonstrate how regularity proper- 
ties of scaling functions can be used to increase the level 
sampling N over the entire x-window during the simu- 
lation. In particular, this will result in an efficient hy- 
brid numerical algorithm that combines the rapidity of 
the FFT for calculating the time evolution, and the inter- 
polation properties of the scaling functions for changing 
the sampling rate. 

The idea is the following. At the beginning of each 
time step, one performs a one level wavelet transform 
of the signal u(t, x) £ Vn from the sampling level N 
to level N — 1, i.e. u(t,x) = ^2 k VN-i,k<PN-i,k(x) + 
y^ /k WN-i,kipN-i,k(x), in order to detect the presence of 
strong gradients. When the maximal absolute value of the 
wavelet coefficients wn-i is greater that a predetermined 
threshold (fixed to 0.005 in our simulations), the num- 



ber of collocation points is doubled using an interpolation 
process. This is done by performing an inverse wavelet 
transform, from level N to level N + 1, on the coefficients 
vn- To this aim, the 2 wavelet coefficients wn are set to 
zero and the result of the inverse transformation is a set 
of 2 N+1 scaling coefficients vn+i that constitute the new 
signal sampling (after multiplying by a \/2 normalization 
factor). 

In the above scheme, an important point to be careful 
about is that no local asymmetry is introduced during the 
resampling process. This is where symmetry properties 
of the scaling functions presented in Section 2 are use- 
ful. Unfortunately, a recomposition using these complex 
scaling functions, slightly couples the real and imaginary 
parts of the signal u(t, x) (for instance, by introducing 
imaginary components in vn+i when the original signal 
is real). This can be avoided, by resampling the signal, us- 
ing a linear combination of the inverse wavelet transform 
and its complex conjugate; in other word, by recomposing 
the signal using the real part of the coefficients au only. 
In fact, we have found this recomposition scheme very 
accurate when done with higher-order wavelets ( J > 8). 
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pling scheme described above [9]. In order to conserve 
the numerical stability, we have also chosen to decrease 
the time step by a factor of 4 after each resampling. The 
maximum wavelet amplitude for the resampling threshold 
has been fixed to 0.005. 

Figure 2 shows a first simulation that leads to strong 
gradients. This is a typical breaking wave phenomena 
that appears in nonlinear optical fibers in the normal dis- 
persion regime, that is for the NTS"— " equation with 
A < 0. In this simulation, A = —900, the initial signal 
is m(0, x) = sech(x) and has been sampled 128 times, the 
initial time range is 7r/32 and the initial time step has 
been set to 7r/1600. The amplitude of the wavelet coeffi- 
cients WN-i(t) is quite enlightning the sharp signal struc- 
ture evolution. Figure 3 shows the evolution of the maxi- 
mal amplitude of these wavelet coefficients as t increases. 
Since each peak is the result of a new sampling, we end 
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Figure 2: Field and wavelet coefficient amplitude for 
a typical breaking wave simulation 

We have written a C-program that simulates (1.1) us- 
ing the split-step Fourier method and the global resam- 



Figure 3: Variation of the maximum amplitude 
wavelet coefficient for the simulation of Figure 2 

up with a final number of 2048 collocation points. A com- 
parison with a non-adaptive simulation (with 2048 sam- 
ples and the corresponding At = 7r/409600) has shown 
that the accuracy of the adaptive simulation is quite sat- 
isfactory. In fact, no difference can be detected on the 
amplitude and phase from the plots. 

The above adaptive algorithm can also be modified to 
allow a decreasing in the number of sampling points when 
sharp gradients tend to disappear. Similarly as above, the 
new signal sampling can be obtained from the scaling co- 
efficients of a one level real- valued wavelet decomposition 
(after multiplying by the normalisation factor l/i/2). 

Interesting solutions that can be used to test the 
"desampling" process are the N-soliton solutions of the 
NLS"+" equation. These are solutions for which the am- 
plitude evolves periodically in time, with a period of 7r/2, 
through successive amplitude peakings. They correspond 
to the initial conditions u(0,x) = N sech(x), where N is 
the solution order (the number of solitons involved in the 
solution). 

Figure 4 shows the variation of the maximal amplitude 



of the wavelet coefficients for the case N = 3, over one pe- 
riod (t = 7r/2) with 128 initial samples and At = 7r/2560. 
The upper threshold has been fixed to 0.005 and the lower 




Figure 4: Variation of the maximum amplitude 
wavelet coefficient for the 3-soliton solution 

one to 30% of the maximal wavelet amplitude after the 
first resampling. This is to avoid a premature decreasing 
of the sampling due to possible oscillatory structures dur- 
ing a resampling relaxation time. The number of sampling 
points varies as 128-256-512-1024-512-1024-512-256-128, 
following the periodicity in the signal amplitude. For this 
particular simulation, we have measured a 10% precision 
on the final signal phase, with respect to the theoretical 
value of 7r/4. On the other hand, the discrepancy of the 
signal amplitude, with respect to the theoretical value, 
turns out to be very small. This is very satisfactory, if we 
take into account that 4 resamplings and 4 "desamplings" 
have been performed during the simulation. 

4. Conclusion and future works 

Wavelets have various fields of applications that go 
well beyond usual signal analysis [10]. Numerical simula- 
tion on wavelet bases is one of them. Here, we have pro- 
posed and studied a simple way to dynamically change 
the sampling rate of a signal that evolves through high 
gradient phenomena. The method was based on the in- 
terpolation properties of highly regular complex symmet- 
ric scaling functions. We have tested our scheme on the 
physically relevant nonlinear Schrodinger model for vari- 
ous high gradient simulations: optical breaking wave, col- 
lapse (not reported here) and bound solitons. In all these 
simulations, the resampling process turned out to be sta- 
ble and very accurate. 

An extension of the present work that is actually un- 
der investigation [9] involves the wavelet-based adaptive 
numerical method developped in Refs. [1, 2]. Briefly, 
the idea is to perform the linear evolution on the wavelet 
spaces rather than the Fourier one, by decomposing the 
signal on . . . Wn-3 ®Wn-2 ®Wn-i • The main advantage 
is that wavelet coefficients that have amplitudes below a 



certain threshold (for instance, 0.0005) can be discarded 
in the numerical algorithm, without introducing signifi- 
cant numerical error. Calculating the evolution on these 
coefficients only, amounts to use a spatial grid that is 
locally adapted to the signal gradients. Our plan is to 
calculate the signal evolution over one time step using 
the same "split-step Fourier algorithm" except that the 
Fourier transform is replaced by a wavelet transform, to- 
gether with the appropriate representation of the second 
derivative operator. The advantage is that nonlinear ef- 
fects can be calculated in the space Vn using a simple 
collocation technique. In addition, we want to keep the 
possibility of globally resampling the signal from Vn to 
Vn+i, using the method described in Section 3, when the 
maximum absolute value of wn-i increases above a fixed 
threshold (for instance, 0.005). 
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